{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "htW5SiGzeXYm"
      },
      "source": [
        "##### Copyright 2018 The TensorFlow Authors.\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "9HGeUNoteaSm"
      },
      "outputs": [],
      "source": [
        "#@title Licensed under the Apache License, Version 2.0 (the \"License\"); { display-mode: \"form\" }\n",
        "# you may not use this file except in compliance with the License.\n",
        "# You may obtain a copy of the License at\n",
        "#\n",
        "# https://www.apache.org/licenses/LICENSE-2.0\n",
        "#\n",
        "# Unless required by applicable law or agreed to in writing, software\n",
        "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "# See the License for the specific language governing permissions and\n",
        "# limitations under the License."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "JJ3UDciDVcB5"
      },
      "source": [
        "# Bayesian Gaussian Mixture Model and Hamiltonian MCMC\n",
        "\n",
        "\u003ctable class=\"tfo-notebook-buttons\" align=\"left\"\u003e\n",
        "  \u003ctd\u003e\n",
        "    \u003ca target=\"_blank\" href=\"https://colab.research.google.com/github/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Bayesian_Gaussian_Mixture_Model.ipynb\"\u003e\u003cimg src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" /\u003eRun in Google Colab\u003c/a\u003e\n",
        "  \u003c/td\u003e\n",
        "  \u003ctd\u003e\n",
        "    \u003ca target=\"_blank\" href=\"https://github.com/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Bayesian_Gaussian_Mixture_Model.ipynb\"\u003e\u003cimg src=\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\" /\u003eView source on GitHub\u003c/a\u003e\n",
        "  \u003c/td\u003e\n",
        "\u003c/table\u003e\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "lin40yCC6eBo"
      },
      "source": [
        "\n",
        "In this colab we'll explore sampling from the posterior of a Bayesian Gaussian Mixture Model (BGMM) using only TensorFlow Probability primitives."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "eZs1ShikNBK2"
      },
      "source": [
        "## Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "7JjokKMbk2hJ"
      },
      "source": [
        "For $k\\in\\{1,\\ldots, K\\}$ mixture components each of dimension $D$, we'd like to model $i\\in\\{1,\\ldots,N\\}$ iid samples using the following Bayesian Gaussian Mixture Model:\n",
        "\n",
        "$$\\begin{align*}\n",
        "\\theta \u0026\\sim \\text{Dirichlet}(\\text{concentration}=\\alpha_0)\\\\\n",
        "\\mu_k \u0026\\sim \\text{Normal}(\\text{loc}=\\mu_{0k}, \\text{scale}=I_D)\\\\\n",
        "T_k \u0026\\sim \\text{Wishart}(\\text{df}=5, \\text{scale}=I_D)\\\\\n",
        "Z_i \u0026\\sim \\text{Categorical}(\\text{probs}=\\theta)\\\\\n",
        "Y_i \u0026\\sim \\text{Normal}(\\text{loc}=\\mu_{z_i}, \\text{scale}=T_{z_i}^{-1/2})\\\\\n",
        "\\end{align*}$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "iySRABi0qZnQ"
      },
      "source": [
        "Note, the `scale` arguments all have `cholesky` semantics. We use this convention because it is that of TF Distributions (which itself uses this convention in part because it is computationally advantageous)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Y6X_Beihwzyi"
      },
      "source": [
        "Our goal is to generate samples from the posterior:\n",
        "\n",
        "$$p\\left(\\theta, \\{\\mu_k, T_k\\}_{k=1}^K \\Big| \\{y_i\\}_{i=1}^N, \\alpha_0, \\{\\mu_{ok}\\}_{k=1}^K\\right)$$\n",
        "\n",
        "Notice that $\\{Z_i\\}_{i=1}^N$ is not present--we're interested in only those random variables which don't scale with $N$.  (And luckily there's a TF distribution which handles marginalizing out $Z_i$.)\n",
        "\n",
        "It is not possible to directly sample from this distribution owing to a computationally intractable normalization term.\n",
        "\n",
        "[Metropolis-Hastings algorithms](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) are technique for for sampling from intractable-to-normalize distributions.\n",
        "\n",
        "TensorFlow Probability offers a number of MCMC options, including several based on Metropolis-Hastings. In this notebook, we'll use [Hamiltonian Monte Carlo](https://en.wikipedia.org/wiki/Hamiltonian_Monte_Carlo)  (`tfp.mcmc.HamiltonianMonteCarlo`). HMC is often a good choice because it can converge rapidly, samples the state space jointly (as opposed to coordinatewise), and leverages one of TF's virtues: automatic differentiation. That said, sampling from a BGMM posterior might actually be better done by other approaches, e.g., [Gibb's sampling](https://en.wikipedia.org/wiki/Gibbs_sampling)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "uswTWdgNu46j"
      },
      "outputs": [],
      "source": [
        "%matplotlib inline\n",
        "\n",
        "from __future__ import absolute_import\n",
        "from __future__ import division\n",
        "from __future__ import print_function\n",
        "\n",
        "import functools\n",
        "\n",
        "import matplotlib.pyplot as plt; plt.style.use('ggplot')\n",
        "import numpy as np\n",
        "import seaborn as sns; sns.set_context('notebook')\n",
        "\n",
        "import tensorflow as tf\n",
        "import tensorflow_probability as tfp\n",
        "\n",
        "tfd = tfp.distributions\n",
        "tfb = tfp.bijectors"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "ovNsKD-OEUzR"
      },
      "outputs": [],
      "source": [
        "def session_options(enable_gpu_ram_resizing=True):\n",
        "  \"\"\"Convenience function which sets common `tf.Session` options.\"\"\"\n",
        "  config = tf.ConfigProto()\n",
        "  config.log_device_placement = True\n",
        "  if enable_gpu_ram_resizing:\n",
        "    # `allow_growth=True` makes it possible to connect multiple colabs to your\n",
        "    # GPU. Otherwise the colab malloc's all GPU ram.\n",
        "    config.gpu_options.allow_growth = True\n",
        "  return config\n",
        "\n",
        "def reset_sess(config=None):\n",
        "  \"\"\"Convenience function to create the TF graph and session, or reset them.\"\"\"\n",
        "  if config is None:\n",
        "    config = session_options()\n",
        "  tf.reset_default_graph()\n",
        "  global sess\n",
        "  try:\n",
        "    sess.close()\n",
        "  except:\n",
        "    pass\n",
        "  sess = tf.InteractiveSession(config=config)\n",
        "\n",
        "reset_sess()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Uj9uHZN2yUqz"
      },
      "source": [
        "Before actually building the model, we'll need to define a new type of distribution.  From the model specification above, its clear we're parameterizing the MVN with an inverse covariance matrix, i.e.,  [precision matrix](https://en.wikipedia.org/wiki/Precision_(statistics%29).  To accomplish this in TF,  we'll need to roll out our `Bijector`.  This `Bijector` will use the forward transformation:\n",
        "\n",
        "- `Y =`  [`tf.matrix_triangular_solve`](https://www.tensorflow.org/api_docs/python/tf/matrix_triangular_solve)`(tf.matrix_transpose(chol_precision_tril), X, adjoint=True) + loc`.\n",
        "\n",
        "And the `log_prob` calculation is just the inverse, i.e.:\n",
        "\n",
        "- `X =` [`tf.matmul`](https://www.tensorflow.org/api_docs/python/tf/matmul)`(chol_precision_tril, X - loc, adjoint_a=True)`.\n",
        "\n",
        "Since all we need for HMC is `log_prob`, this means we avoid ever calling `tf.matrix_triangular_solve` (as would be the case for `tfd.MultivariateNormalTriL`). This is advantageous since `tf.matmul` is usually faster owing to better cache locality.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "nc4yy6vW-lC_"
      },
      "outputs": [],
      "source": [
        "class MVNCholPrecisionTriL(tfd.TransformedDistribution):\n",
        "  \"\"\"MVN from loc and (Cholesky) precision matrix.\"\"\"\n",
        "\n",
        "  def __init__(self, loc, chol_precision_tril, name=None):\n",
        "    super(MVNCholPrecisionTriL, self).__init__(\n",
        "        distribution=tfd.Independent(tfd.Normal(tf.zeros_like(loc),\n",
        "                                                scale=tf.ones_like(loc)),\n",
        "                                     reinterpreted_batch_ndims=1),\n",
        "        bijector=tfb.Chain([\n",
        "            tfb.Affine(shift=loc),\n",
        "            tfb.Invert(tfb.Affine(scale_tril=chol_precision_tril,\n",
        "                                  adjoint=True)),\n",
        "        ]),\n",
        "        name=name)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "JDOkWhDQg4ZG"
      },
      "source": [
        "The `tfd.Independent` distribution turns independent draws of one distribution, into a multivariate distribution with statistically independent coordinates. In terms of computing `log_prob`, this \"meta-distribution\" manifests as a simple sum over the event dimension(s).\n",
        "\n",
        "Also notice that we took the `adjoint` (\"transpose\") of the scale matrix. This is because if precision is inverse covariance, i.e., $P=C^{-1}$ and if $C=AA^\\top$, then $P=BB^{\\top}$ where $B=A^{-\\top}$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Pfkc8cmhh2Qz"
      },
      "source": [
        "Since this distribution is kind of tricky, let's quickly verify that our `MVNCholPrecisionTriL` works as we think it should."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 153
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 507,
          "status": "ok",
          "timestamp": 1538760458123,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "GhqbjwlIh1Vn",
        "outputId": "d9c48cc8-5dd2-45c9-9ccc-3be0da3619ba"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "true mean: [ 1. -1.]\n",
            "sample mean: [ 1.00028062 -1.00010502]\n",
            "true cov:\n",
            " [[ 1.0625   -0.03125 ]\n",
            " [-0.03125   0.015625]]\n",
            "sample cov:\n",
            " [[ 1.06412733 -0.03126175]\n",
            " [-0.03126175  0.01559312]]\n"
          ]
        }
      ],
      "source": [
        "def compute_sample_stats(d, seed=42, n=int(1e6)):\n",
        "  x = d.sample(n, seed=seed)\n",
        "  sample_mean = tf.reduce_mean(x, axis=0, keepdims=True)\n",
        "  s = x - sample_mean\n",
        "  sample_cov = tf.matmul(s, s, adjoint_a=True) / tf.cast(n, s.dtype)\n",
        "  sample_scale = tf.cholesky(sample_cov)\n",
        "  sample_mean = sample_mean[0]\n",
        "  return [\n",
        "      sample_mean,\n",
        "      sample_cov,\n",
        "      sample_scale,\n",
        "  ]\n",
        "\n",
        "dtype = np.float32\n",
        "true_loc = np.array([1., -1.], dtype=dtype)\n",
        "true_chol_precision = np.array([[1., 0.],\n",
        "                                [2., 8.]],\n",
        "                               dtype=dtype)\n",
        "true_precision = np.matmul(true_chol_precision, true_chol_precision.T)\n",
        "true_cov = np.linalg.inv(true_precision)\n",
        "\n",
        "d = MVNCholPrecisionTriL(\n",
        "    loc=true_loc,\n",
        "    chol_precision_tril=true_chol_precision)\n",
        "\n",
        "[\n",
        "    sample_mean_,\n",
        "    sample_cov_,\n",
        "    sample_scale_,\n",
        "] = sess.run(compute_sample_stats(d))\n",
        "\n",
        "print('true mean:', true_loc)\n",
        "print('sample mean:', sample_mean_)\n",
        "print('true cov:\\n', true_cov)\n",
        "print('sample cov:\\n', sample_cov_)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "N60z8scN1v6E"
      },
      "source": [
        "Since the sample mean and covariance are close to the true mean and covariance, it seems like the distribution is correctly implemented. Now, we'll use `MVNCholPrecisionTriL` and  stock`tfp.distributions` to specify the BGMM prior random variables:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "xhzxySDjL2-S"
      },
      "outputs": [],
      "source": [
        "dtype = np.float32\n",
        "dims = 2\n",
        "components = 3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "xAOmHhZ7LzDQ"
      },
      "outputs": [],
      "source": [
        "rv_mix_probs = tfd.Dirichlet(\n",
        "    concentration=np.ones(components, dtype) / 10.,\n",
        "    name='rv_mix_probs')\n",
        "\n",
        "rv_loc = tfd.Independent(\n",
        "    tfd.Normal(\n",
        "        loc=np.stack([\n",
        "            -np.ones(dims, dtype),\n",
        "            np.zeros(dims, dtype),\n",
        "            np.ones(dims, dtype),\n",
        "        ]),\n",
        "        scale=tf.ones([components, dims], dtype)),\n",
        "    reinterpreted_batch_ndims=1,\n",
        "    name='rv_loc')\n",
        "\n",
        "rv_precision = tfd.Wishart(\n",
        "    df=5,\n",
        "    scale_tril=np.stack([np.eye(dims, dtype=dtype)]*components),\n",
        "    input_output_cholesky=True,\n",
        "    name='rv_precision')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 88
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 244,
          "status": "ok",
          "timestamp": 1538760469863,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "KSTp8aAIAv0O",
        "outputId": "9b18cf98-8514-471b-a56f-429be5b30d3d"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "tfp.distributions.Dirichlet(\"rv_mix_probs/\", batch_shape=(), event_shape=(3,), dtype=float32)\n",
            "tfp.distributions.Independent(\"rv_loc/\", batch_shape=(3,), event_shape=(2,), dtype=float32)\n",
            "tfp.distributions.Wishart(\"rv_precision/\", batch_shape=(3,), event_shape=(2, 2), dtype=float32)\n"
          ]
        }
      ],
      "source": [
        "print(rv_mix_probs)\n",
        "print(rv_loc)\n",
        "print(rv_precision)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "8ZOG0OR815Nr"
      },
      "source": [
        "Using the three random variables defined above, we can now specify the joint log probability function. To do this we'll use `tfd.MixtureSameFamily` to automatically integrate out the categorical $\\{Z_i\\}_{i=1}^N$ draws."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "CpLnRJr2TXYD"
      },
      "outputs": [],
      "source": [
        "def joint_log_prob(observations, mix_probs, loc, chol_precision):\n",
        "  \"\"\"BGMM with priors: loc=Normal, precision=Inverse-Wishart, mix=Dirichlet.\n",
        "\n",
        "  Args:\n",
        "    observations: `[n, d]`-shaped `Tensor` representing Bayesian Gaussian\n",
        "      Mixture model draws. Each sample is a length-`d` vector.\n",
        "    mix_probs: `[K]`-shaped `Tensor` representing random draw from\n",
        "      `SoftmaxInverse(Dirichlet)` prior.\n",
        "    loc: `[K, d]`-shaped `Tensor` representing the location parameter of the\n",
        "      `K` components.\n",
        "    chol_precision: `[K, d, d]`-shaped `Tensor` representing `K` lower\n",
        "      triangular `cholesky(Precision)` matrices, each being sampled from\n",
        "      a Wishart distribution.\n",
        "\n",
        "  Returns:\n",
        "    log_prob: `Tensor` representing joint log-density over all inputs.\n",
        "  \"\"\"\n",
        "  rv_observations = tfd.MixtureSameFamily(\n",
        "      mixture_distribution=tfd.Categorical(probs=mix_probs),\n",
        "      components_distribution=MVNCholPrecisionTriL(\n",
        "          loc=loc,\n",
        "          chol_precision_tril=chol_precision))\n",
        "  log_prob_parts = [\n",
        "      rv_observations.log_prob(observations), # Sum over samples.\n",
        "      rv_mix_probs.log_prob(mix_probs)[..., tf.newaxis],\n",
        "      rv_loc.log_prob(loc),                   # Sum over components.\n",
        "      rv_precision.log_prob(chol_precision),  # Sum over components.\n",
        "  ]\n",
        "  sum_log_prob = tf.reduce_sum(tf.concat(log_prob_parts, axis=-1), axis=-1)\n",
        "  # Note: for easy debugging, uncomment the following:\n",
        "  # sum_log_prob = tf.Print(sum_log_prob, log_prob_parts)\n",
        "  return sum_log_prob"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "QM1idLJazkGC"
      },
      "source": [
        "Notice that this function internally defines a new random variable. This is necessary since the `observations` RV depends on samples from the RVs defined further above."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "7jTMXdymV1QJ"
      },
      "source": [
        "## Generate \"Training\" Data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "rl4brz3G3pS7"
      },
      "source": [
        "For this demo, we'll sample some random data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "1AJZAtwXV8RQ"
      },
      "outputs": [],
      "source": [
        "num_samples = 1000\n",
        "true_loc = np.array([[-2, -2],\n",
        "                     [0, 0],\n",
        "                     [2, 2]], dtype)\n",
        "random = np.random.RandomState(seed=42)\n",
        "\n",
        "true_hidden_component = random.randint(0, components, num_samples)\n",
        "observations = (true_loc[true_hidden_component] +\n",
        "                random.randn(num_samples, dims).astype(dtype))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "zVOvMh7MV37A"
      },
      "source": [
        "## Bayesian Inference using HMC"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "cdN3iKFT32Jp"
      },
      "source": [
        "Now that we've used TFD to specify our model and obtained some observed data, we have all the necessary pieces to run HMC.\n",
        "\n",
        "To do this, we'll use a [partial application](https://en.wikipedia.org/wiki/Partial_application) to \"pin down\" the things we don't want to sample. In this case that means we need only pin down `observations`. (The hyper-parameters are already baked in to the prior distributions and not part of the `joint_log_prob` function signature.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "tVoaDFSf7L_j"
      },
      "outputs": [],
      "source": [
        "unnormalized_posterior_log_prob = functools.partial(joint_log_prob, observations)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "a0OMIWIYeMmQ"
      },
      "outputs": [],
      "source": [
        "initial_state = [\n",
        "    tf.fill([components],\n",
        "            value=np.array(1. / components, dtype),\n",
        "            name='mix_probs'),\n",
        "    tf.constant(np.array([[-2, -2],\n",
        "                          [0, 0],\n",
        "                          [2, 2]], dtype),\n",
        "                name='loc'),\n",
        "    tf.eye(dims, batch_shape=[components], dtype=dtype, name='chol_precision'),\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "TVpiT3LLyfcO"
      },
      "source": [
        "### Unconstrained Representation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "JS8XOsxiyiBV"
      },
      "source": [
        "Hamiltonian Monte Carlo (HMC) requires the target log-probability function be differentiable with respect to its arguments.  Furthermore, HMC can exhibit dramatically higher statistical efficiency if the state-space is unconstrained.\n",
        "\n",
        "This means we'll have to work out two main issues when sampling from the BGMM posterior:\n",
        "\n",
        "1. $\\theta$ represents a discrete probability vector, i.e., must be such that $\\sum_{k=1}^K \\theta_k = 1$ and $\\theta_k\u003e0$.\n",
        "2. $T_k$ represents an inverse covariance matrix, i.e., must be such that $T_k \\succ 0$, i.e., is [positive definite](https://en.wikipedia.org/wiki/Positive-definite_matrix).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Vt9SXJzO0Cks"
      },
      "source": [
        "To address this requirement we'll need to:\n",
        "\n",
        "1. transform the constrained variables to an unconstrained space\n",
        "2. run the MCMC in unconstrained space\n",
        "3. transform the unconstrained variables back to the constrained space.\n",
        "\n",
        "As with `MVNCholPrecisionTriL`, we'll use [`Bijector`s](https://www.tensorflow.org/api_docs/python/tf/distributions/bijectors/Bijector) to transform random variables to unconstrained space.\n",
        "\n",
        "- The [`Dirichlet`](https://en.wikipedia.org/wiki/Dirichlet_distribution) is transformed to unconstrained space via the [softmax function](https://en.wikipedia.org/wiki/Softmax_function).\n",
        "\n",
        "- Our precision random variable is a distribution over postive semidefinite matrices. To unconstrain these we'll use the `FillTriangular` and `TransformDiagonal` bijectors.  These convert vectors to lower-triangular matrices and ensure the diagonal is positive. The former is useful because it enables sampling only $d(d+1)/2$ floats rather than $d^2$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "_atEQrDR7JvG"
      },
      "outputs": [],
      "source": [
        "unconstraining_bijectors = [\n",
        "    tfb.SoftmaxCentered(),\n",
        "    tfb.Identity(),\n",
        "    tfb.Chain([\n",
        "        tfb.TransformDiagonal(tfb.Softplus()),\n",
        "        tfb.FillTriangular(),\n",
        "    ])]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "0zq6QJJ-NSPJ"
      },
      "outputs": [],
      "source": [
        "[mix_probs, loc, chol_precision], kernel_results = tfp.mcmc.sample_chain(\n",
        "    num_results=2000,\n",
        "    num_burnin_steps=500,\n",
        "    current_state=initial_state,\n",
        "    kernel=tfp.mcmc.TransformedTransitionKernel(\n",
        "        inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(\n",
        "            target_log_prob_fn=unnormalized_posterior_log_prob,\n",
        "            step_size=0.065,\n",
        "            num_leapfrog_steps=5),\n",
        "        bijector=unconstraining_bijectors))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "_ceX1A3-ZFiN"
      },
      "outputs": [],
      "source": [
        "acceptance_rate = tf.reduce_mean(tf.to_float(kernel_results.inner_results.is_accepted))\n",
        "mean_mix_probs = tf.reduce_mean(mix_probs, axis=0)\n",
        "mean_loc = tf.reduce_mean(loc, axis=0)\n",
        "mean_chol_precision = tf.reduce_mean(chol_precision, axis=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "kmpTFZcVmByb"
      },
      "source": [
        "Note: through trial-and-error we've predetermined the `step_size` and `num_leapfrog_steps` to approximately achieve an [asymptotically optimal rate of 0.651](https://arxiv.org/abs/1001.4460). For a technique to do this automatically, see the examples section in `help(tfp.mcmc.HamiltonianMonteCarlo)`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "QLEz96mg6fpZ"
      },
      "source": [
        "We'll now execute the chain and print the posterior means."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "3B2yJWVmNcrm"
      },
      "outputs": [],
      "source": [
        "[\n",
        "    acceptance_rate_,\n",
        "    mean_mix_probs_,\n",
        "    mean_loc_,\n",
        "    mean_chol_precision_,\n",
        "    mix_probs_,\n",
        "    loc_,\n",
        "    chol_precision_,\n",
        "] = sess.run([\n",
        "    acceptance_rate,\n",
        "    mean_mix_probs,\n",
        "    mean_loc,\n",
        "    mean_chol_precision,\n",
        "    mix_probs,\n",
        "    loc,\n",
        "    chol_precision,\n",
        "])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 306
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 257,
          "status": "ok",
          "timestamp": 1538760907401,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "bqJ6RSJxegC6",
        "outputId": "25a86519-852f-4368-f8d4-2914ff1db494"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "    acceptance_rate: 0.666\n",
            "      avg mix probs: [ 0.3831611   0.25915343  0.3576861 ]\n",
            "\n",
            "            avg loc:\n",
            " [[-1.87364268 -1.79061651]\n",
            " [-0.02801551  0.0025092 ]\n",
            " [ 1.88933122  1.92743731]]\n",
            "\n",
            "avg chol(precision):\n",
            " [[[ 1.01128054  0.        ]\n",
            "  [-0.06019538  0.97167695]]\n",
            "\n",
            " [[ 1.2244916   0.        ]\n",
            "  [ 0.21678646  1.02939916]]\n",
            "\n",
            " [[ 0.98610836  0.        ]\n",
            "  [-0.11749489  0.96770507]]]\n"
          ]
        }
      ],
      "source": [
        "print('    acceptance_rate:', acceptance_rate_)\n",
        "print('      avg mix probs:', mean_mix_probs_)\n",
        "print('\\n            avg loc:\\n', mean_loc_)\n",
        "print('\\navg chol(precision):\\n', mean_chol_precision_)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 286
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 11557,
          "status": "ok",
          "timestamp": 1538760926241,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "zFOU0j9kPdUy",
        "outputId": "187b368d-0fb7-4e49-88b5-e215f4a27d9d"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAENCAYAAAASUO4dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGQhJREFUeJzt3X9s1PXhx/HXlVKBzrXXAoVrZYBFm+KAYltYKiStWzZZ\nBBIHosucbv7BD4U5UDMXfjgY6AZm/Ai4uAkmuGkwysz4Y1sMfjMF5Je4CFQGFlxb2kKPVmTQo+37\n+we76931rj+Pfu59fT4SEu7X5/PuJ8eTd9/3uTuXMcYIAGCtJKcHAADoHUIOAJYj5ABgOUIOAJYj\n5ABgOUIOAJYj5EhYf/rTn1RSUqLJkyersbEx5Laqqirl5eWptbX1po6hrKxM+/fvv6n7AAg5YiY8\nWnv27FFxcbEOHz4cCOfkyZM1efJk3XPPPZo/f7727dvXbhsTJ07U5MmTVVBQoMmTJ2vNmjXdHktz\nc7NefPFFbd++XUePHlVaWlq7+7hcru7/kEAcSnZ6AEhM77zzjl588UW98sormjhxoqqqquRyuXTk\nyBG5XC7V19drz549WrRokVauXKnZs2cHHvv73/9eU6dO7dX+L168KJ/Pp9tvv723P8pN09LSogED\nBjg9DCQAZuSIuTfffFO/+c1v9Oqrr2rixIkht/nfSJyZmalHHnlETz75pH77299GvE9nfD6ffv3r\nX2vatGmaPn261q5dq+vXr+vs2bO67777JElFRUV69NFHO91WXV2dFixYoClTpui73/2udu3aFbit\ntbVVL7/8sr7zne/o7rvv1gMPPKDa2tqI29m9e7fKyso0depUvfzyyyG3bdmyRYsXL9bTTz+twsJC\nvfPOO/rXv/6lefPmqaioSNOmTdPq1avV3NwsSdq8eXPgt5Hm5mYVFBRo/fr1kqSmpiZNmDBBly9f\nls/n09NPP60pU6aoqKhIc+bMkdfr7dIxRIIwQIyUlpaaJ5980pSUlJjPPvss5LbKykqTl5dnWlpa\nQq7/4osvzJ133mnOnDkT2Ma+ffu6tL/f/e535sEHHzRer9d4vV7z4IMPmo0bN4bsr7W1NeJjw8fz\nwx/+0PzqV78yPp/PnDx50kydOtXs37/fGGPMK6+8Yu6//35z9uxZY4wx5eXlpqGhod02//3vf5tJ\nkyaZw4cPG5/PZ9atW2fGjx8f+Hk2b95sxo8fb9577z1jjDFNTU3m+PHj5pNPPjGtra2mqqrKzJgx\nw7z22mvGGGP2799v7r//fmOMMUePHjXf/va3zdy5c40xxuzbt8/MmjXLGGPMG2+8YebPn2+amppM\na2urOX78uPnqq6+6dAyRGJiRI6b27duniRMn6o477ujS/bOysiQp5MXIRYsWqbi4WEVFRSouLg6Z\nHQf761//qkWLFsntdsvtduuJJ57Q7t27JbXN6k0XZvfnz5/Xxx9/rGXLlmngwIHKy8vTnDlz9Je/\n/EWS9NZbb+mpp57SN77xDUnSnXfeGXHN/W9/+5vKysp09913a+DAgVqyZEm7+xQUFKisrEySlJKS\novz8fE2YMEEul0sej0dz587VoUOHAvc9d+6cGhsbdejQIf3gBz9QbW2trl69qsOHD6uoqEiSlJyc\nrIaGBlVUVMjlcik/P1+pqamd/txIHKyRI6aef/55bd26Vc8995zWrl3b6f39SxTp6emB67Zu3dql\nNfK6ujp5PJ7AZY/HowsXLkjq3guZFy5cUFpamgYPHhyyrePHj0uSampqdNttt3VpPCNGjAhcHjx4\ncMjPJSnkdkk6e/asXnjhBX366ae6du2aWlpaNH78eEnSLbfcorvuuksHDx7U4cOHtWDBApWXl+vI\nkSM6ePCgHnnkEUnSrFmzVFNTo5///Oe6fPmyZs6cqaeeeor1936EGTliKiMjQzt27NCRI0e0atWq\nTu//97//XUOHDtWYMWMC13VlFi3dmM1XVVUFLldXV2v48OHdHvPw4cPV2Nio//73v4Hrzp8/H9jW\niBEj9MUXX3S6nWHDhqmmpiZw+erVq2poaAi5T/h/MKtWrdLYsWP1j3/8Q4cPH9bPfvazkJ+/sLBQ\nBw4c0MmTJ/XNb35ThYWF+uCDD/Tpp5+qsLBQ0o0Z+aJFi7Rnzx698cYb2rt3b+A3E/QPhBwxN2zY\nML322mv64IMPtG7dusD1xphApOrr67Vz505t3bpVS5cu7dF+ZsyYoW3btsnr9crr9Wrr1q2aNWtW\nyP464r99xIgRKigo0EsvvSSfz6fy8nK99dZbmjlzpiRpzpw52rhxo86dOydJ+uyzz9qdly5J3/ve\n97R3714dPXpU169f16ZNmzr9Ga5cuaKvfe1rGjx4sM6cOaM///nPIbcXFxdr9+7dys3NVXJysqZM\nmaJdu3YpJydHbrdbkvTRRx/p1KlTam1t1ZAhQ5ScnMxsvJ9haQUxEzzbHDFihHbs2KEf/ehHGjRo\nkObOnSuXy6WioiIZYzRkyBDddddd2rRpk0pKSkK2s2DBAiUltc0xSkpKtHnz5nb7W7hwoa5cuaKZ\nM2fK5XLpvvvu0/z58yOOp7PxbtiwQStXrtS0adOUlpamJUuW6Fvf+pYk6bHHHtP169f1k5/8RA0N\nDRo7dqy2bNnSbp08NzdXK1as0NKlS3X16lU99thjgdcAonn22We1fPly/eEPf1B+fr6+//3v68CB\nA4HbCwoK1NTUFFgPz83N1aBBgwKXpRunWq5cuVK1tbVKTU3VjBkzAv8JoX9wma7+HgsAiEssrQCA\n5Qg5AFiOkAOA5Qg5AFiOkAOA5Rw7/bC6utqR/Xo8Hsf2HW84Fm04FqF6ejxSh3b/DVnourSUyMlm\nRg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5vrMT\nACJo8Pm6df/0lJSbNJLOEXIA+J/weF9saox4v6G3pLW7LvixfR11Qg6g3wuOcHC8yxtqQ+6Xl54V\ncp9IQfdvry9jTsgB9EvRZt/h8T5xqVKSlO/OUXlDbSDm8aTXIa+vr9eWLVvU0NCgpKQk3XvvvZox\nY0YsxgYAMdXZ0ok/4v54Bwu+LnhmHm1W3pd6HfIBAwboxz/+sUaPHq1r167p2Wef1cSJE5WdnR2L\n8QFATPgj3tHSSXCsz1yqCLntdveYkMd1NDO3bo08PT1d6enpkqRBgwYpOztbXq+XkAOIC5HWvzua\neYcHPPj6291jdOJSpfLdOVH358TZKzFdI6+rq9O5c+c0bty4WG4WALotWsD98Y4WbL/K+lOSpJzM\nO27SCGPHZYwxsdjQtWvXtGrVKj3wwAMqKiqKxSYBWKbR1+z0ECS1X0YJD7g/0uf/Ux54zMjb8qJu\nzx9z//JKvjsnsLQSvEZ+s2fj0b58OSYz8paWFm3YsEHTp0/vcsSd+sZyvi29DceiDcciVE+PR+rQ\n4TdhNN0TLeL+gPvjXVVbH3hMdlZm4PqOgi5FjriTbwaSYvQW/W3btiknJ4ezVQDElUgRr6qtD4m4\nFBr14Fm6FDobD14bj5eISzGYkZeXl+uf//ynRo0apWeeeUYul0sPPfSQJk2aFIvxAUCXRTozJZqm\nsw26ZXR6h/eJFPG89Ky4irgUg5Dn5eXpzTffjMVYAKDHOlpSiSRRIi7xzk4ACSD8jT7h54f7X9yM\nJjsrM/D3kbfltYt4vK2JhyPkAKwWaSYezH+miT/mwdEOFingkuI+4hIhB2CpjgIe6Y0+HZ0P3pWA\nS/EZcYmQA0gAHb1TU7oR6DOXKqLGvKPzw+M13sEIOQBrBZ+dEi3ifsGflRLOpmWUSAg5AOt099t7\nJEX9fJTgD7+yYRklEkIOICHku3OizsqDIx7+qYW2xjsYIQdgraG3pOliU2OnX/YQbdYt2RvvYIQc\ngNX8Ye4o6IkY72CEHEBC6Mo39SRawP0IOYCElqjxDkbIASSk/hBwP0IOwDrpKSntTkHsT+EOR8gB\nWKk/hztcTL5YAgDgHEIOAJYj5ABgOUIOAJYj5ABgOUIOAJYj5ABgOUIOAJYj5ABgOUIOAJYj5ABg\nOT5rBUDMXLlY5/QQYs7j8ai6utrpYUiS0jyeiNczIwcAyxFyALAcIQcAyxFyALAcIQcAyxFyALAc\nIQcAyxFyALAcIQcAyxFyALAcIQcAyxFyALBcTD40a9u2bTp69KjS0tK0fv36WGwSANBFMZmRl5aW\n6pe//GUsNgUA6KaYhDwvL0+pqamx2BQAoJtYIwcAyxFyALCcY98Q5InyTReJvu94w7Fow7EIxfFo\nE+/HImYhN8bIGNPl+zv11Unx9LVNTuNYtOFYhOJ4tImnYxHtP5SYhHzjxo06ceKELl++rAULFmju\n3LkqLS2NxaYBAJ2ISciXLFkSi80AAHqAFzsBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAs59g7OxEf\nvjYsy+khxIUvr7c4eiy+ulDr2L5hP2bkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4A\nliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AliPkAGA5Qg4AluPL\nl4EovvQ1t7vu6yn8k0H8YUaOfu1LX3O7YEe6Lvg2IN4wvUC/FCnekXh9jcpISeuLIQE9RsjR7wRH\n2+trlKR2sfZfH+k+X/qaWWJBXOHZiH7FH/HgUEe6LEkVl+sDfx9za2bI7JyYI57wTES/ER7x4FB3\npuJyvcbcmhlxe36EHU7hmYeEF2kppeJyvU431nZ5G7lpWYHHs2aOeEPIkbCiBVySTjfW6rOGyk63\ncWd6TuDv/lk5MUe8IeRISJ3Nwv0Rr26siLoNT9qYkNjnpmURc8QlQo6EEx7x4Fm41D7iNRc+D3n8\niGFjQ24PDrp/iSXaflknhxN41iGhRHtBM3gpJTzgdTVnI24rOOjhs/PwFz4BJ8Uk5MeOHdOOHTtk\njFFpaalmz54di80C3RIp4h3NwqMFvK7mrIaPGK2aC58HYt6dMTArR1/r9TOutbVVf/zjH7VixQq5\n3W794he/UFFRkbKzs2MxPqBLokW8o4BX1XrbbSc7K6MPRgvEVq8/a+X06dMaOXKkhg0bpuTkZJWU\nlOjQoUOxGBvQJeHnc4efHx4e8apar6pqvWr9z+XAH79Icfe7Mz1HuWlZHS6rMBuHE3odcq/Xq8zM\ntid2RkaGvN7o/xiAWIp0doqkdueIB0dcUki8I12OtuwSzddTkok4HHNTnnkul+tmbBYIEent9sGz\n8c8aKlXdWNEu4t3lSRsTmI1HQsDhtF4/AzMyMnTx4sXAZa/XK7fb3enjPB5Pb3fdY07uGzdP+GmG\nHUm67dYObx8+YnTE6/3LKrE+hzwen5PxOCanxPux6HXIc3NzVVNTowsXLsjtduvDDz/UkiVLOn1c\ndXV1b3fdIx6Px7F9x5t4f3L2RvhsPBL/C5vBM/XgFzv9Z6x0ZW28t+LtOcm/kzbxdCyi/ZvtdciT\nkpL005/+VGvWrJExRmVlZcrJyen8gUAfi3ZGSvD14TNxT9qYdvcPno2zrIJ4EJNn4aRJk7Rx48ZY\nbArosq+nJOtLX7MyUtLk9TVqzK2ZET/RcPiI0YFZebQlE79I540Hz8aJOOIRz0QkBH/M/YI/7EqK\nHGip7d2d/tuDZ+Dh2yDiiFc8G5FQ/DPn04217UIcSfjSSfBj/GepjLk1k4gjrvGMhLU6+iLk3LSs\nQMzDP662o8AHx9uPiCPe8ayElTqKePhaeXDMwyMe6dxwIg7b8MxEwghfJw8WaRYeHvHw0wuJOGzB\nsxMJJTjm/uWVYNHenRltFi4RccQ/nqGwkv/Uw2j8yyuRZt3hpygScdiOZymsFRzZ4KgHn1ceSbTr\niThsxTMVCSF8ht7ZZ6EEr6UTcNiOZywSRmfLLcGihZ6Iw0Y8a5FQCDH6o15/sQQAwFmEHAAsR8gB\nwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKE\nHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAsR8gBwHKEHAAs\nl9ybBx84cEC7du1SZWWl1q1bp7Fjx8ZqXOgjX12odXoIccHj8ai6utrpYQA90qsZ+ahRo7Rs2TLl\n5+fHajwAgG7q1Yzc4/HEahwAgB5ijRwALNfpjHz16tVqbGwMXDbGyOVyad68eSosLOzxjp2czfOb\nRBuORRuORSiOR5t4Pxadhnz58uU3ZcdOvbDEi1ptOBZtOBahOB5t4ulYRPsPhaUVALBcr0J+8OBB\nLViwQKdOndILL7ygtWvXxmpcAIAu6tVZK8XFxSouLo7VWAAAPcDSCgBYjpADgOUIOQBYjpADgOUI\nOQBYjpADgOUIOQBYjpADgOUIOQBYjpADgOUIOQBYjpADgOV69aFZvZE6dLgj+230NTu2785cuVjn\n9BAAWIgZOQBYjpADgOUIOQBYjpADgOUIOQBYjpADgOUIOQBYzrHzyG3W4POFXE5PSXFoJADAjLxb\nGny+dhH3Xw8ATiHkXdRZrKNFHgBuNpZWOtCTMAc/hiUXAH2BkIcJj/fFpsaI9xt6S1pfDAcAOkXI\nFT3e5Q21Ee+fl54VuE+0oDMbB9BX+n3I/REPj/eJS5UR75/vzlF5Q63y0rNCHufHTB1AX+u3IY8U\n8BOXKnXmUoUq60+F3Dcn846Qy+ExBwAn9buQBy+jXGxqDJmB/9/pv+n8f8qjPjYn8w6duVQhiZgD\niB/95vTD8NMDwyN+5lKFzv+nXFW19SF/JAXiHj5TD+dfVmF9HEBf6hcz8vCAS21LKZICM3F/uJvO\nNkiSbhmd3uV9EHEATukXIfeLth4eKeLBRt6WJ+nG0srt7jHKd+dIUmBZhYgDcFLCh7yzFzWDIy6F\nzsKzszIDf/dH3C8vPSvkDBUiDsApCRvyjpZTgiMeTXZWZshM3C/fndPuBU4iDsBJCRfySAGXIr+5\nZ+RteRFj7o94cMCDl1QkllMAxI+ECnn4Mko0/iWSyvpT/5t1h8bcH3H//VgTBxDPEirkXZHvztGJ\nS5W63T1Gt7vH6MylinZv+Oks4BIRBxA/EibknX1SYbQ37vjD7v97pMcwAwcQz3oV8p07d+rIkSNK\nTk5WVlaWFi5cqCFDhsRqbD0W7fNO/Esu4VGPFHlm3wBs0auQT5gwQQ8//LCSkpL0+uuva/fu3Xr4\n4YdjNbaYixT44PX08NsJOAAb9DrkfuPGjdNHH33U6wH1VHh0u/qlEJHiTsAB2CRmn7Wyd+9eFRQU\nxGpzvdZZjNNTUkLu479MxAHYxmWMMR3dYfXq1WpsbFt+MMbI5XJp3rx5KiwslCS9/fbb+vzzz7Vs\n2bIu77jR19zDIXddtFl5vMY6LSVhXnsG0Ic6DXln3n//fb333ntasWKFBg4c2OXH9UXI/Wz5Hs0r\nF+v6dH8ej0fV1dV9us94xbEIxfFoE0/HwuPxRLy+V1PAY8eO6d1339Xzzz/frYj3tXiONwD0Vq9C\n/uqrr6q5uVlr1qyRdOMFz8cffzwmAwMAdE2vQr5p06ZYjQMA0EP95huCACBREXIAsBwhBwDLEXIA\nsBwhBwDLEXIAsBwhBwDL9fot+gAAZzEjBwDLEXIAsBwhBwDLEXIAsBwhBwDLEXIAsFy/+26xnTt3\n6siRI0pOTlZWVpYWLlyoIUOGOD0sxxw4cEC7du1SZWWl1q1bp7Fjxzo9pD537Ngx7dixQ8YYlZaW\navbs2U4PyTHbtm3T0aNHlZaWpvXr1zs9HEfV19dry5YtamhoUFJSku69917NmDHD6WFFZvqZTz75\nxLS0tBhjjNm5c6d5/fXXHR6Rs6qqqkx1dbVZtWqVOXPmjNPD6XMtLS3miSeeMHV1deb69etm2bJl\nprKy0ulhOebkyZOmoqLCLF261OmhOO7SpUumoqLCGGPM1atXzeLFi+P2udHvllYmTJigpKQbP/a4\nceNUX1/v8Iic5fF4NHLkSKeH4ZjTp09r5MiRGjZsmJKTk1VSUqJDhw45PSzH5OXlKTU11elhxIX0\n9HSNHj1akjRo0CBlZ2fL6/U6O6go+l3Ig+3du1cFBQVODwMO8nq9yszMDFzOyMiI23+scE5dXZ3O\nnTuncePGOT2UiBJyjXz16tVqbGwMXDbGyOVyad68eSosLJQkvf322xowYIDuuecep4bZZ7pyPNDG\n5XI5PQTEkWvXrumll17So48+qkGDBjk9nIgSMuTLly/v8Pb3339fH3/8sVasWNFHI3JWZ8ejP8vI\nyNDFixcDl71er9xut4MjQjxpaWnRhg0bNH36dBUVFTk9nKj63dLKsWPH9O677+qZZ57RwIEDnR4O\nHJabm6uamhpduHBBzc3N+vDDD/v9bynGGBk+S0/SjbN4cnJy4vdslf/pd59+uHjxYjU3N+vWW2+V\ndOMFz8cff9zhUTnn4MGD2r59u7788kulpqZq9OjReu6555weVp86duyYtm/fLmOMysrK+vXphxs3\nbtSJEyd0+fJlpaWlae7cuSotLXV6WI4oLy/XypUrNWrUKLlcLrlcLj300EOaNGmS00Nrp9+FHAAS\nTb9bWgGAREPIAcByhBwALEfIAcByhBwALEfIAcByhBwALEfIAcBy/w/9HHUX7Gu42QAAAABJRU5E\nrkJggg==\n",
            "text/plain": [
              "\u003cmatplotlib.figure.Figure at 0x7f40ea6065d0\u003e"
            ]
          },
          "metadata": {
            "tags": []
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "ax = sns.kdeplot(loc_[:,0,0], loc_[:,0,1], shade=True)\n",
        "ax = sns.kdeplot(loc_[:,1,0], loc_[:,1,1], shade=True)\n",
        "ax = sns.kdeplot(loc_[:,2,0], loc_[:,2,1], shade=True)\n",
        "plt.title('KDE of loc draws');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "NmfNIM1c6mwc"
      },
      "source": [
        "## Conclusion"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "t8LeIeMn6ot4"
      },
      "source": [
        "This simple colab demonstrated how TensorFlow Probability primitives can be used to build hierarchical Bayesian mixture models."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "Bayesian_Gaussian_Mixture_Model.ipynb",
      "provenance": [],
      "toc_visible": true,
      "version": "0.3.2"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
